arXiv:1503.07973vl [stat.ME] 27 Mar 2015 


Accelerated least squares estimation for systems 
of ordinary differential equations 

Itai Dattner Shota Gugushvili 

March 30, 2015 


Abstract 

We study the problem of parameter estimation for a system of ordinary dif¬ 
ferential equations based on noisy observations on the solution of the system. A 
classical estimation approach to this problem is the least squares method. Owing to 
a highly nonlinear character of the least squares criterion function and the need to 
employ repetitive numerical integration of the system, the latter method becomes 
computationally intense for most realistic systems. We propose the accelerated 
(ACCEL) least squares method, which starts from a preliminary y^-consistent es¬ 
timator of the parameter of interest and next through a Newton-Raphson type step 
turns it into an asymptotically equivalent estimator to the least squares estimator. 
Additional computational burden of this step is negligible. We demonstrate excel¬ 
lent practical performance of the ACCEL least squares approach via simulations 
and real data examples. The method enables the researcher to obtain both point 
and interval estimates. 

March 30, 2015 


1 Introduction 


Systems of ordinary differential equations (ODEs) are commonly used for the math¬ 
ematical m odeling of the rate of cha nge of dynamic processes (e.g., in mathematical 
biology, s ee Edelstein-Keshet ( 2005h : in the theory of chem ical reaction networks, see 
Eeinbergl ( 1979 ) and Sontagl l 200 ih : and in biochemistry, see Voil ( 2000h ). The systems 
we have in mind take the form 


\x'{t)=F{x{t),d), fG[0,l], 

I -*^(0) = ^ 


( 1 ) 


where x{t) = {x\{t) ^. ,Xii[t)y is a c/-dimensional state variable, 0 = (0i,...,0p)' 
denotes a p-dimensional parameter, while the column c/-vector x(0) = ^ defines the 
initial condition. We define rj := (^,9) and denote the solution to ([T]) corresponding to 
the parameter rj by 

x{ri,t) := {xi{ri,t),...,xd{ri,t)f. 
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Knowledge regarding the system parameters ^ and 6 is of vital importance for 
the study of a process that ([T]i models. Indeed, these parameters affect the qualita¬ 
tive properties of the system, and their knowledge allows one to predict the system 
behaviour. However, in practice the parameter 0, and possibly also are unknown 
to the researcher. Typically they cannot be measured directly, but have to be inferred 
from noisy measurements of the process under study. Let rjo = (^o, ^o) be the ‘true’ 
parameter value that governs the underlying process. The common statistical model 
considered for the noisy measurements of the process at time instances f i,..., f„ (deter¬ 
ministic or random; not necessarily equally spaced), is the additive measurement error 
model, 

Yij=Xi{Tlo,tj) + £ij, i=l,...,d,j=l,...,n, ( 2 ) 

where the random variables e,j are independent measurement errors (not necessarily 
Gaussian). Based on observation pairs i = \ the goal is to 

estimate the parameter tjq. 

A classical approach to parameter estimation for ordinary differential equations is 
the nonlinear least squares (NLS) method. Its use is based on the observation that the 
problem at hand in its essence is a nonlinear regression problem, where the regression 
function x(t 7 ,-) is defined implicitly as the solution to ([T]i. The least squares estimator 
t?n = {^n,Sn) of Tjo is defined as a minimiser of the least squares criterion function 


Tj„ = e„) = argmin^ £ £ (7,^ 

;=i;=i 

=: argmin^ Rniri). 


(3) 


The strongest justification for the use of the least s quar es e s timato r lies in its attrac¬ 
tive asymptotic properties; see e.g.. |jennrich ( 1969 ) and (1981). In most practical 
applications the solution x{rj,-) to ([T]) is nonlinear in the parameter 77, and therefore 
some iterative procedure has to be used to compute the nonlinear least squares estima¬ 
tor. Such procedures require an initial guess for a minimiser rjn and then proceed by 
constructing successive approximations to the least squares estimator (in a direction 
guided by the gradient of the criterion function, when a gradient based optimisation 
method, e.g. the Levenberg-Marquardt method, is used). However, the noisy and non¬ 
linear character of the optimization problem may lead for the procedure to end up in 
a local minimum of the least squares criterion function, especially when good initial 
guesses of the parameter values are not available. Furthermore, in most applications 
the system (ID is nonlinear and does not have a closed form solution. In that case at 
every step of the iterative procedure one has to numerically integrate (ID (as well as 
the system of the associated sensitivity equations in order to compute the gradient of 
the criterion function, in case a gradient-based optimisation method is used). Since 
the number of iterations made until convergence of the algorithm can be ascertained is 
usually large, in most cases this leads to a computational bottleneck. This is the case 
especially in mathematical biology and biochemistry, where a highly nonlinear charac¬ 
ter of dependence of the solution x( 77,-) on the parameter 77 leads to ‘stiff’ integra tion 
problems. For a penetrat ing discussion of these point see e.g. iRamsav et al.l (l2007h and 


Voit and Almeidal (2004). 
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A Bayesian approach to parameter estimation fo r ordinary diff e rentia l equations 
is a vi a ble alternative to the NLS metho d; see, e.g., Gelman et al. ( 19961) . Girolami 
( 2008h . Girolami and Calderhead ( 2011 ). However, it will not be discussed in the 
present work. We only remark that like the NLS method, its application also requires a 
considerable computational effort. 

In order to byp ass computational problem s ass ociated with c omputation of the least 
squares estimator. iBellman and RothI ( 1971 ) and Varah ( 1982h suggested to estimate 
the solution x(t7o, •) u sing nonparametric techni q ues (e .g., splines). This approach was 
stu died rigorously in Gugushvili and Klaassen (l2012h (other relevant re ferences are 


tea rigorously in Kiugusnyili ana Klaassen i QZUl^) (otner relevant re terences are 
Brunell (l2008l) . IVuiacic et alJ (l2014l) and Dattner and Klaass^ ( 201-5h ') and works 


e-g- 

as follows: the observations are first smoothed, which results in an estimator x„(-) 
for the solution ^(tjo, •) of the system, and by differentiation, in an estimator x^(-) for 
y (tjo, •) ■ Then the estimator for 9q is the minimizer over 6 of the smooth function 


/' 




(4) 


where w is an appropriate weight function, and || • || denotes the standard Euclidean 
norm. Hence, this approach bypasses the need to integrate the system numerically, and 
as a result the parameter estimates can be computed e xtremely quickly, especially when 
E in ([T]i is linear in 9. Under regularity conditions iGugushvili and KlaassenI (12012h 
show that this smooth and match estimator (SME) 9„ has the of convergence 

to 0. By the general statistical theory, the ^^-rate of convergence is in fact the best 
rate one can expect in the present context. This result thus puts the smooth and match 
method on a solid theoretical ground. 

Note that execution of this method does not require the knowledge of the initial 
values. However, it cannot be used to estimate initial values. If estimation of initial 
values is of interest, then once the estimator 0„ is at hand, one may obtain an estimator 
by minimizing with respect to ^ the criterion 


[ ||x„(f)-^- [ F{xn{s)-9„)ds\\^dt. 

Jo Jo 


Notice that this is a linear least squares optimization problem and hence is easy to 
execute. 

The SME described above is criticized for not being statistically efficient. In infor¬ 
mal terms this means that the resulting estimators do not squeeze as much information 
out of the data as the least squares estimator does (in more formal terms, their asymp¬ 
totic variance is larger than that of the least squares es timator). Hence, sometimes it 
is suggested (see, e.g.. lSwartz and BremermannI (Il975h for an early reference) to use 
this method only for generating preliminary estimates that should be used later as an 
initial guess for more accurate methods, that, unfortunately, are not computationally 
efficient. In this paper we present a new algorithm, ACCEL (accelerated) least squares, 
that enjoys both the computational efficiency of the SME approach and the statisti¬ 
cal efficiency of the least squares estimators. The method is studied theoretically in 
an unpublished report by Gugushvili and Klaassen, while here we will concentrate on 
numerical aspects. 
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The rest of the paper is organized as follows: in Section|2]we describe the ACCEL 
least squares algorithm. Section [3] presents a simulation study illustrating the perfor¬ 
mance of our method, while Section |4] contains numerical results based on real data 
examples. Sectionj^summarizes our contribution and outlines potential future research 
directions. 


2 Accelerated least squares estimator 


When one adpots an asymptotic point of view on statistics, all the estimators with the 
same asymptotic variance can be considered as equivalent. We now demonstrate how 
once some ^^-consistent estimator rj„ of the parameter Tj is available (e.g. the SME), 
one can obtain an asymptotically equivalent estimator to the least squares estimator 
in just one extra easy step. The method i s referred to as the one-step method in the 
statistical literature, see e.g. Section 5.7 in ivan der Vaartl(ll998h . 

Consider the estimating equations 




(5) 


j=i 


where 

Wriit^y) = {y-x{ri,t)){x^iri,t)f, 

with jcjj (tj,?) denoting the derivative of x{rj,t) with respect to Tj. Up to omission of a 
constant multiplier —2, these estimating equations are obtained by differentiating the 
least squares criterion function R„ (tj ) from (O with respect to tj . A value of tj = (^, 0) 
that solves (|5]) is called a Z-estimator. In general there can be multiple solutions to (|5]l, 
one of them being the least squares estimator. Since tj is a p -f -dimensional vector, 
we will need p + d estimating equations, so that ‘Th(tj) =0 stands for p + d equations. 

The accelerated (ACCEL) least squares estimator Tj„ of tjq is defined as a solution 
in TJ of the equation 

- Vn) = 0 . 

This corresponds to replacing 'Th(tj) by its tangent line at Tj„, and ij„ is then the ‘tj- 
intercept’ of the tangent line. If is invertible, the estimator ij„ can be ex¬ 

pressed as 

In = ln- 


If Tjn is ‘close’ to the solution of the estimating equations (|5]l, the ACCEL estimator tj,, 
will also be ‘close’ to it, and hopefully, even ‘closer’ than rj„. In the numerical analysis, 
when looking for a zero of an equation, such a procedure is iterated several times, each 
time taking the previously found value as a new start ing point, and is known u nder 
the name of Newton’s method, see e.g.. Section 2.3 in Burden and EairesI (2001). On 
the other hand, the asymptotic theory of statistics tells us that given certain regularity 
assumptions on the statistical model, the estimator ij„ is just as accurate asymptotically 
as a Z-estimator obtained from solving 'Tn(T7) = 0. Thus from the asymptotic point of 
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view there is no difference between the least squares estimator and the estimator Tj,,; 
both are equally accurate. At the same time, the extra computational cost needed to 
evaluate the ACCEL estimator rj^ is negligible. Indeed, the only essential additional 
thing that has to be done for its computation in comparison to computation of the 
SME is the evaluation of and the derivative matrix This reduces 

to requiring just one numerical integration of the sensitivity and variational equations 
associated with O, as we will now illustrate. 

In what follows, it is helpful to think of E in ([1]) as a function of rj rather than 
only 0. Thus, we write the right-hand side E of ([Til as F{x{rj,t),rj). Differentiating 
both sides of ([T]i with respect to 77 and interchanging the order of a f-derivative with an 
T]-derivative, we get 

This system is a matrix differential equations and is usually referred to in the literature 
as a system of sensitivity equations. By replacing 77 with fj„ we arrive at the system 

V( 0 ) = (l, 0 )^ 

where we have defined s(t) := -^xirinj). Observe that x(rj„,-) is a known function, 

because it can be found by integrating ([T]i for parameter values and 9„. Consequently, 
the system of sensitivity equations is a linear system with time-dependent coefficients, 
and hence is relatively straightforward to integrate. 

By differentiating (jb]) one more time with respect to 77 and replacing 77 with rj„ we 
arrive at the following set of variational equations: 


{ 1 - 2(0 = F^'ri (4Vn,t),T]„))+ F”^{x{^„ ,t),r]„)s{t) 

+ J7«)J?«>(0} ^(0 

+Kix{rin,t),rin)z{t), 

2 ( 0 ) = 0 , 


where we have set z{t) := -^x(r],t). Eor each z,-, i = l,...,c/, the system (O is a 
matrix differential equation and again is a linear system with time-varying coefficients. 
Here also we can treat x and s as known, for they can be obtained through numerical 
integration of ([T]i and (|7]i. The process of obtaining variational equations can be made 
automatic through a software implementation, and we provide such a Matlab code. 

Integration of ([T]i, (0 and (|8]l for the parameter value 77 ,, will allow us to compute 
'i'nifjn) and and consequently, the ACCEL estimator 77 ,,. Note that nu¬ 

merical integration of the variational equations (or at least the sensitivity equations) is 
usually required when computing the least squares estimator via gradient-based opti¬ 
mization meth ods (unless the gradient is availa ble analytically) or using the Bayesian 
approach (see iGirolami and CalderheadI (1201 ih l. However, we need to do this only 
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Remark 1 A seemingly more general nonautonomous system than the autonomous 
system o, 

r 7{t) = F{x{t),t-e), t e [0,1], 

I ^( 0 ) = ^, 

may and will be reduced to o by a simple substitution x{t) = (x^(f),f)^, t G [0,1], 
and ^ = (^^,0)^. 


3 Simulation study 


3.1 Linear ODE 


We start with illustrating the performance of the ACCEL estimator when used to esti¬ 
mate the parameter and initial value of a one-dimensional linear ordinary differential 
equation 

ix'{t) = dQx{t), t G [o,r], 

U(0) = ^o- 


(9) 


This is a toy example, but it allows us to explore the practical performance of the 
ACCEL method in great detail and to compare it to the theoretically expected results. 
Advanced examples will be considered later on. 

Given observations T,/s, the ACCEL method requires first to have at hand a ^/n- 
rate estimator of 9q and ^o- As mentioned in the previous sections, the SME provides 
us with such an estimator. However, this method is based on estimating the derivative 
x', which is hard to do accurately in practice. In the case where the symbol of the 
system of ODEs is linear in functions of the parameter 0, as is the case in (|9]l, one can 
avoid the estimation of derivatives and use an integral S ME. Indeed, in such cases one 
can use some version of the so called ’ integr al approach’ (IHimmelblau et al.l (1196711 1 as 
was studied in iDattner and KlaassenI (l2015h . Eor systems linear in the parameter 0 as 
in (|9]l the idea works as follows. 

Note that in (|9]l F(x(t);9) = g(x(t))9 holds, where the measurable function g : 

—>■ maps the c/-dimensional column vector x into a dx p matrix. Let x„{-) be 

an estimator of x(rio,-), and denote G„(t) = fQg(x„(s),s)ds, A „ = Jq G„(t)dt, B„ = 
Tif GI (t)G„(t)dt, and be the d x d identity matrix. Then Dattner andKlaassen 
( 20151) show that the direct estimators 


t = ^ (/,,-A„B-iG^(f))x„(f)dr, (10) 

On = Gl{t)(^Ut)-l,)dt, (11) 


are ^^-consistent. In case the initial value is known, (fTTT) may be used with 
replaced by Besides the required statistical properties, the extensive simulation 
study presented in the aforementioned paper suggests that this approach is much more 
accurate in finite samples comparing to the derivative SME resulting from minimizing 
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(IHi. Thus, we use the integral SME (fTOl i- (fTTl i whenever applicable, and the derivative 
SME (|4]i otherwise. 

We choose to estimate the solution x using local polynomial estimators, which are 
consistent and ‘automatically’ correct for the boundaries. Under the assumption that 
X are C“-functions for som e real a > 1, we will approximate them by polynomials of 
degree i = [aj as follows ( Tsvbako^2009l Section 1.6). Let 


Uiu) = (i,m,mV(2!),...,m 7(^0)^> weM, 

v(f) = ^x{t),x'{t)b,x"{t)b^,...,x^^\t)b^'^ , f € K., 


where b = b„ > 0 is a bandwidth, the {£ -b 1)-vector U{u) is a column vector, and v(f) 
isad X {(. + l)-matrix. Let /f(-) be some appropriate kernel function and define 


v„(f) = arg min 






The local polynomial estimator of order £ of x{t) is the first column of the dx (.^ -b 1)- 
matrix v„(f), i.e.,x„(f) = v„(f)t/(0). 

In the sequel we apply the estimation procedure described above to a set of band- 
widths B ;= {/^min, • • • ,bmax}, and for a given b G B we denote the resulting ACCEL 
parameter estimator by fj„ h. We then choose 


— Bn,b^ 


( 12 ) 


where 


^ = argminj^ L (5); 

dGB 


i=\j=\ 


Here Xi{fj„ h,-) is the solution of (ffl) using the ACCEL estim ator 

The theoretical analysis of iDattner and KlaassenI (12015h shows that in order to have 
the Y^-rate for the parameter estimators, once should take a bandwidth b = 

Thus, we set B = x (ci, ...,Cjv), N a positive integer, and the c/s depend on the 
grid of points on which we evaluate the kernel estimator. Last, we use local estimators 
polyno mials of order 1, with K{t) = 3/4(1 — f^)l{|f| < 1} (cf., Dattner and Klaass^ 
( 2015l) l. where !{•} stands for the indicator function. Other kernels are also possible. 

The solution of the initial value problem (|9]l is x(f) =xoexp(0oO- We generate 
(pseudo) random observations from the model 


T, = ^oexp(0of,7) + eA 

where tj G {0(0.1)10} {n = 101), and ej ^ A(0,0.05^), j = 1... ,n. We consider 
00 e {-1,-0.8,-0.6,-0.4,-0.2,0.1,0.3,0.5,0.7,0.9} 
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and ^0 G {0.5,1}. For each pair {^o, 9o) we run a Monte Carlo study of 500 samples 
of Fi... ,Fioi, where in each sample we apply the ACCEL least squares method as 
described above, and the nonlinear least squares. We note that we observed some 
setups, where the ACCEL estimator may perform differently, depending whether initial 
values are estimated or not. Thus, in this example we used scaling of the time t, such 
that t € [0,1], which stabilized the procedure. However, scaling was not needed in 
all other examples we consider below. This simulation study enables us to estimate 
the asymptotic variance of the least squares and the ACCEL least squares methods. 
We then compare the results to the true asymptotic variance. The true and estimated 
asymptotic variances can be obtained for each set of parameters and initial values by 
inverting the Eisher information matrix. Assuming the observation model (|2l), where 
we have homoscedastic errors, and assuming that the observed time points 
are uniformly distributed over [0,7’], the Eisher information matrix is given by the 
{d + p)x {d + p) matrix 

where we can obtain s{t) := -^x{r\j) as in (|2l). However, the Eisher information 
matrix depends on the true values of the parameters, initial values, and which are 
not known in practice. So in general, we will construct a fully data driven confidence 
interval by estimating the Eisher information matrix. We estimate <7^ by 


^2 = 


1 


d{n — 1) 




1=1./=i 


where stands for the solution of the system (|9]) using the estimated parame¬ 

ters and initial values using the ACCEL method. Then an estimate for the asymptotic 
variance of the parameter rjj is given by /yy^(Tj„)/«,where stands for the y'th 

diagonal element of the inverse Eisher information matrix evaluated in point fj,,. When 
i(-) has no closed form, the integral in (fOl ) will be evaluated using numerical integra¬ 
tion (in the subsequent examples we use the trapezoidal rule). 

A direct computation gives that in model (l9|l the asymptotic variance of depends 
on 9, but is independent of the values of ^ itself. In Eigure[T]we see the estimated vari¬ 
ance of the ACCEL estimators (plus signs) and that of the NLS (circles), for estimating 
^0 based on 500 simulation runs. The estimates are superimposed on the theoretical 
asymptotic variance (dashed line). The top plot is for = 0.5 and the bottom one is 
for 1*0 = 1. As the theory suggests, independently of the values of ^, the true asymp¬ 
totic variance is the same. Note that in this specific numerical example the estimated 
variances of the ACCEL and NLS estimators are the same. This is not surprising, since 
in order to apply the NLS we used as the initial point in the parameter space the SME 
(resulted from using a kernel with the largest bandwidth; this choice was arbitrary). 
The estimated variances agree with the asymptotic one. We note that the grid of 9o 
does not include 0, where the asymptotic variance equals zero. 

In Eigure|2]we see similar plots corresponding to estimating the asymptotic vari¬ 
ances of 9„. Here the variance has different order, depending on the value of ^o- Again, 
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Figure 1: The estimated variance of the ACCEL (plus signs) and NLS (circles) es¬ 
timators and respectively, based on 500 simulations with n = 101 and ej ~ 
A(0,0.05^), j = 1... ,n. The estimates are superimposed on the theoretical asymptotic 
variance (dashed line). The upper plot is for i^o = 0.5 and the lower one is for = 1. 

the estimated variances of the ACCEL (plus signs) and NLS (circles) estimators are 
the same, and both agree with the asymptotic one (dashed line). Similar plots were 
obtained when considering other values for ( 7 ^, and therefore we do not present them 
here. 

The discussion above implies that we can provide confidence intervals together 
with the point estimates. An approximate I — a level confidence interval for rjj is 
given by 


[f}j,„ - Zl-alj/ {flj,n)/n, fin - Zl-alj/ (i?j,n)/«], (14) 

where Zi-a is the 1 — a quantile of the standard normal distribution. In Table [T] we 
present the coverage of this confidence interval based on a Monte Carlo study with 
500 simulations for different experimental setups. The results should be compared to 
the nominal coverage of 95%. We consider 4 setups denoted by A,B,C,D according 
to (^0 = 1/2 ,00 =-1), (^0 = 1/2, 00 = 1 ), (^0 = 1 ,00 =-1), (^0 = 1 ,00 =-1), re¬ 
spectively. Each scenario is tested for n =21, and n = 51. Table [T] presents the point 
and interval estimates for the parameters of each scenario. We see that the coverage 
of the data driven confidence intervals is satisfying across the different experimental 
scenarios. 

We close this example by the following illustration. In Table |2] we present for 
both ACCEL and NLS methods, the square root of the average of the estimates of the 
asymptotic variance, over the 500 simulations (denoted by ’ASYM’). Next to that we 
present standard errors of the point estimates as calculated based on the 500 simulations 
(denoted by ’STL’). The results agree with each other. 


3.2 Lotka-Volterra system 

The Lotka-Volterra system of ODEs dEdelstein-Keshel (2005)) is a population dynam¬ 
ics model that describes evolution over time of the populations of two species, predators 
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Figure 2: The estimated variance of the ACCEL (plus signs) and NLS (circles) es¬ 
timators On and On, respectively, based on 500 simulations with ii = 101 and ej ^ 
N {0,0.05^), j = I... ,n. The estimates are superimposed on the theoretical asymptotic 
variance (dashed line). The upper plot is for = 0.5 and the lower one is for = 1. 


Table 1: Means of point estimates and actual coverage of interval estimates for 
the parameters of model (|9l) according to 4 different experimental setups. The re¬ 
sults are based on 500 simulation runs. The observations are generated according to 
Yj = ^oexp{Ootj) + Ej, where tj G {0(0.5)10} {n = 21 ), or tj G {0(0.2)10} (n = 51 ) 
and Sj ^ A(0,0.05^), j = The point estimates are given by (fT^ : the interval 

estimates are defined in (fT4l) 


Setup 




ACCEL 

Mean Coverage 

Mean 

NLS 

Coverage 

n=21 

A 

^0 

0.500 

0.501 

0.942 

0.501 

0.942 



00 

-1.000 

-1.002 

0.946 

-1.002 

0.946 


B 

^0 

0.500 

0.500 

0.928 

0.500 

0.928 



00 

1.000 

1.000 

0.938 

1.000 

0.938 


C 

^0 

1.000 

0.999 

0.932 

0.999 

0.932 



00 

-1.000 

-0.997 

0.940 

-0.997 

0.940 


D 

^0 

1.000 

1.000 

0.944 

1.000 

0.944 



00 

1.000 

1.000 

0.948 

1.000 

0.948 

n=51 

A 

^0 

0.500 

0.500 

0.944 

0.500 

0.944 



00 

-1.000 

-0.998 

0.944 

-0.998 

0.944 


B 

^0 

0.500 

0.500 

0.946 

0.500 

0.946 



00 

1.000 

0.999 

0.958 

0.999 

0.958 


C 

^0 

1.000 

0.999 

0.932 

0.999 

0.932 



00 

-1.000 

-0.999 

0.938 

-1.000 

0.938 


D 

^0 

1.000 

1.000 

0.948 

1.000 

0.948 



00 

1.000 

1.001 

0.952 

1.001 

0.952 
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Table 2: Standard errors of the point estimates as calculated based on the 500 simula¬ 
tions (denoted by ’STE’). Square root of the average of the estimates of the asymptotic 
variance, over the 500 simulations (denoted by ’ ASYM’). 


Setup 




ACCEL 

STE ASYM 

NLS 

STE ASYM 

n=21 

A 

^0 

0.500 

0.027 

0.031 

0.027 

0.031 



00 

-1.000 

0.121 

0.140 

0.121 

0.140 


B 

^0 

0.500 

0.016 

0.012 

0.016 

0.012 



00 

1.000 

0.045 

0.034 

0.045 

0.034 


C 

^0 

1.000 

0.027 

0.034 

0.027 

0.034 



00 

-1.000 

0.061 

0.079 

0.061 

0.079 


D 

^0 

1.000 

0.016 

0.016 

0.016 

0.016 



00 

1.000 

0.022 

0.023 

0.022 

0.023 

n=51 

A 

^0 

0.500 

0.017 

0.017 

0.017 

0.017 



00 

-1.000 

0.080 

0.076 

0.080 

0.076 


B 

^0 

0.500 

0.010 

0.009 

0.010 

0.009 



00 

1.000 

0.028 

0.026 

0.028 

0.026 


C 

^0 

1.000 

0.017 

0.019 

0.017 

0.019 



00 

-1.000 

0.040 

0.043 

0.040 

0.043 


D 

^0 

1.000 

0.010 

0.010 

0.010 

0.010 



00 

1.000 

0.015 

0.014 

0.015 

0.014 


and their preys. The system takes the form 

f x\(t) = 9ixi(t)-92Xi(t)x2(t), 

\ x' 2 (t) =- 93 x 2 ( 1 )+ 94 x 1 (t)x 2 (t). 

Here xi represents the size of the prey population and X 2 of the predator population. 
In Table [ 3 ] we see the empirical coverage of the 95% confidence intervals based on a 
Monte Carlo study consisting of 500 simulation runs for different sample sizes. The 
experimental setup is as follows; the observed time points are equidistant on [0,10]; the 
errors are normal with zero mean and standard deviation a — 0.05; the initial values are 
^0 = (1,1 /2)^, and the parameters are 9q = (1 /2,1 /2,1 /2,1 /2)^. The point estimates 
are given by (fT2l i. while the interval estimates are defined in (fT4l i. As expected, the 
coverage is much better when the sample size is larger. The performance of the ACCEL 
and NLS methods is similar. 

Last, for both ACCEL and NLS methods, in Table |4| we present the square root 
of the average of the estimates of the asymptotic variance, over the 500 simulations 
(denoted by ’ASYM’). Next to that we present standard errors of the point estimates 
as calculated based on the 500 simulations (denoted by ’STE’). The results agree with 
each other. 
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Table 3; Means of point estimates and actual coverage of interval estimates for the pa¬ 
rameters of model (fTSl l. where the initial values are ^ = (Ij 1 /2)^, and the rate param¬ 
eters are 0o = (1 /2, 1 /2,1 /2, 1 /2)^. The results are based on running 500 simulations. 
The observed time points are equidistant on [0,10], and the errors are normal with zero 
expectation and (7 = 0.05. The ACCEL point estimates are given by (fT2li : the interval 
estimates are defined in (fT4l) . 


Setup 



ACCEL 

Mean Coverage 

Mean 

NTS 

Coverage 

n=21 


1.000 

1.000 

0.932 

0.999 

0.928 


ki 

0.500 

0.500 

0.936 

0.500 

0.934 


01 

0.500 

0.502 

0.942 

0.501 

0.942 


02 

0.500 

0.502 

0.932 

0.501 

0.938 


03 

0.500 

0.500 

0.910 

0.501 

0.916 


04 

0.500 

0.500 

0.918 

0.501 

0.922 

n=51 


1.000 

1.000 

0.958 

1.000 

0.966 


^2 

0.500 

0.500 

0.954 

0.500 

0.948 


01 

0.500 

0.502 

0.964 

0.500 

0.968 


02 

0.500 

0.501 

0.968 

0.500 

0.964 


03 

0.500 

0.500 

0.958 

0.500 

0.958 


04 

0.500 

0.500 

0.952 

0.500 

0.958 


Table 4: Standard errors of the point estimates as calculated based on the 500 simula¬ 
tions (denoted by ’STE’). Square root of the average of the estimates of the asymptotic 
variance, over the 500 simulations (denoted by ’ ASYM’). 




ACCEL 

NTS 

Setup 


STE 

ASYM 

STE 

ASYM 

n=21 

^1 

0.025 

0.023 

0.025 

0.023 



0.020 

0.019 

0.020 

0.019 


01 

0.027 

0.026 

0.027 

0.026 


02 

0.022 

0.021 

0.022 

0.021 


03 

0.022 

0.020 

0.022 

0.020 


04 

0.020 

0.018 

0.020 

0.018 

n=51 

^1 

0.015 

0.016 

0.014 

0.016 


& 

0.013 

0.013 

0.013 

0.013 


01 

0.016 

0.017 

0.016 

0.017 


02 

0.013 

0.014 

0.013 

0.014 


03 

0.013 

0.013 

0.013 

0.014 


04 

0.012 

0.012 

0.012 

0.012 
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4 Real data examples 


In this section we study several realistic examples, also employing real data. To check 
applicability of our method, the emphasis will be on examples with small or moderate 
sample sizes. 


4.1 Nitrogene oxide reaction 

The system 

= 0i(126.2 —x(f))(91.9 —.r(f))^— 92{x{t)Y, 
\x[0)=0 

describes the reversible homogeneous gas phase reaction of nitrogene oxide, 

2 NO + O 2 ^2N02. 


The forward reaction is the third order, while the reverse reaction is the second order. 
Time in (fThl l is measured in minutes and the variable x(t) m odels pressure fall at time f. 
For additional chemic al background see Bodensteinl(ll922h . Based on the experimental 
data from Table 39 in Bodenstei^ (ll92 ^, pa r ameters of equation (fThll were estim ated 


via different methods in 


Bellman et al.l (1967 ): Van Dom selaar and Hemke J ( 1975h . 


see 


pp. 18-19; Esposito and FloudasI ( 200ol) . Section 7.4: iKim and Shend (l2010h . Section 
3.1; Tioa and Biegleil ( 1991 ). Problem 6 on p. 381; and Varahl ( 1V82 ). see pp. 37- 
38. The results obtained in these papers are summarised in Table l5p| We also remark 
that this problem is one of the six test pro blems in param eter estimation for ordinary 
differential equations that were included in FloudasI ( 1999 ). 


Table 5: Parameter estimates for model (fTbl l obtained in the literature. 


Paper 

Estimate of 0i 

Estimate of 02 

Bellman et al. ('19671 

0.4577 X 10^^ 

0.2797 X 10-5 

Van Domselaar and Hemker 119751 

0.45 X 10-5 

0.27 X 10-5 

Esposito and Eloudas 120001 

0.4593 X 10-5 

0.28285 X 10-5 

Kim and Sheng 120101 

0.46 X 10-5 

0.28 X 10-5 

Tioa and Bieeler 119911 

0.4604 X 10-5 

0.2847 X 10-5 

Varah 119821 

0.46 X 10-5 

0.27 X 10-5 


Our interest in this example first went in the following direction: we used the real¬ 
istic estimated parameter values from the literature, generated an artificial set of data 
from (fTbl l and checked how well the ACCEL estimator performs in this case. We 

'Note that IVaraJ ll982h gives five different pai'ameter estimates con'esponding to different values of 
the smoothin g parameter u sed in his method. Of these estimates we report only the first pair and refer 
to Table 4 in ivarahl fl982ll for the remaining ones. Note also that lEsnosito and FloudasI j2Q0Qh use two 
approaches (collocation method and integration method in their terminology) and with the second of them 
identify another local soluti on to the problem, namely 9i — 0.1306 x 10“^,02 = 0.90393 (see Table 11 in 
lEsposito and FloudasI j2000h l. which we did not report in Table[5] 
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Figure 3: The solution to (fTSl l is given by the solid line, while the observations are 
indicated by pluses. The parameter is Tjo = = (0,0.4577 x 10^^,0.2797 x 

10^^)^, and the observations were generated uniformly over tj € {0(2)40}, {n = 21). 
The errors ej were generated from A^(0, (7^) with mean zero and = 0.25. 


also present the estimation results using the nonlinear least squares estimator. Accord¬ 
ingly, w£_took_the_parameter estimates 0i = 0.4577 x 10^^ and 02 = 0.2797 x 10^^ 
from iBellman et al.l (119671) together with the initial condition ^ = x(0) = 0, thus rjo = 
(^,01,02)^. Then we generated observations uniformly over fy € {0(2)40}, (n = 21), 
according to (|2|i, where the i.i.d. measurement errors ej were generated from the nor¬ 
mal distribution A(0, <7^) with mean zero and variance (7^ = 0.25. The solution to (fTbl l 
with these parameter values is plotted in Figure [3] with a solid line, while the observa¬ 
tions are indicated by pluses. This setup was chosen to mimic the real data scenario 
related to this model, as described later on. 

The fact that 0i and 02 are small numbers, combined with the fact that their mag- 
nitudes are rather different, renders their estimation a difficult task, cf. p. 1303 in 
Esposito and Floudas ( 2000l) . In Table |6] we see the empirical average of point esti¬ 
mates and the empirical coverage of interval estimates based on Monte Carlo study con¬ 
sisting of 500 runs. The point estimates are given by (fTSl i. while the int erval estimates 
are de f ined in (fT4ll. We note that when est i mating 0 = (0|, 09), un like iBeUmanet^ 
( 1967 1. Van Domselaar and Hemker ( 1975 ). lTioa and Bieglei ( 1991 ) and IVarahId 1982 ). 
we did not assume that the initial condition x(0) = 0 was known, but estimated it as 
well. Notice also that our method exploits linearity in the parameters an d therefore it is 


not req uired to provide it with an initial guess in the parameter space (in lBellman et al 

and 02 = 10^^ were used). 


( 1967h and other related papers the initial guesses 0i = 10 ® ^^ ’ 


We see that even with a small sample as 21 observations, the point and interval esti¬ 
mates are satisfying, and again, we do not observe a substantial difference between the 
ACCEL and NLS methods. 

We next te s ted ou r approach on the real data for the model (fTbl l given in Table 39 


BodensteinI ( 1922h and reproduced in Table I in iBellman et al. ( 1967 ). There are 


in total 14 observations available on the interval [0,39], excluding the initial condition 
x{0) = 0.E This time we did not estimate the initial condition and considered it to be 


^Note ttiat in Table 39 in iBodensteinI 1 19221) and in Table I in lBellman et all 1 19671) the observation 48.8 
corresponding to the time instance t — \9 appears to contain a typo: we tentatively corrected it to 38.8. The 
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Table 6: Means of point estimates and actual coverage of interval estimates for the 
parameters of model (fTbl l. where the initial value ^ is zero and the parameters are 
01 = 0.4577 X 10^^ and 02 = 0.2797 x 10^^. The results are based on 500 simulation 
runs. There are 21 observations given on a uniform grid on [0,40], and the errors are 
normal with zero expectation and = 0.25. The ACCEL point estimates are given by 
(O; the interval estimates are defined in (fl4li . 


Setup 



ACCEL 

Mean Coverage 

NLS 

Mean Coverage 

n=21 

^1 

0 

1.491e-02 

0.938 

7.960e-03 

0.942 


01 

4.577e-06 

4.576e-06 

0.954 

4.577e-06 

0.952 


02 

2.797e-04 

2.788e-04 

0.932 

2.798e-04 

0.930 


zero, which agrees with the physical phenomenon the model describes. The estimation 
results are displayed in Table|7] Both point and interval estimates obtained through the 
ACCEL and NLS methods are presented. The results are very similar. 


Table 7: P oint estimates for t he parameters of model (fTbl l based on the real data of 
Table 39 in Bodenstehil ( 1922h . We consider the initial value to be zero. The ACCEL 


point estimates are given by (fT2li : the confidence intervals were generated according to 


dl. The left and right interval points are denoted by CI(L) and CI(R), respectively. 



Point 

ACCEL 

CI(L) 

CI(R) 

Point 

NLS 

CI(L) 

CI(R) 

01 

4.579e-06 

4.255e-06 

4.903e-06 

4.577e-06 

4.253e-06 

4.901e-06 

02 

2.791e-04 

1.923e-04 

3.658e-04 

2.796e-04 

1.928e-04 

3.665e-04 


A comparison to the results given in Table|5]shows that this is essentially the same 
result as already reported in the literature using the least squares estimator: this illus¬ 
trates the fact that just one iteration in the Newton-Raphson type procedure suffices to 
arrive at an asymptotically equivalent estimator to the least squares estimator of the pa¬ 
rameter of interest, provided a preliminary estimator is already within the range 
of the true parameter. In Eigure|4|we plot the data from iBellman et alJ (119671) and the 
solution to (flbl l evaluated with ACCEL fitted values of 0i and 02. The fit appears to be 
satisfactory given a simplistic character of the model ( fTSl l. 


4.2 Barnes’ problem 


This p roblem was studied in IVarahl ( 1982 ). who refers to Van Domselaar and Hemkei 


(119751) for exact experimental details. The system under study represents chemical 
reaction equations and is actually a variant of the Lotka-Volterra system that was treated 


same correction was applied in Table 24 in lEsposito and FloudasI l2000l) and in Table 1 in iKim and Shend 
120 id) . 
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Figure 4: The solution to (fThl l is given by the solid line, while the observ ations are indi¬ 


cated b y pluses. The parameters were estimated using the real data from iBellman et al 
( 19671) . The initial value is considered to be known and equals zero. 


above. The system takes the form 


= 0lXl{t) - 92Xi{t)x2{t), 
x'2{t) = e2Xi{t)x2{t) - 92,X2{t). 


(17) 


In the real data of Varah ( 1982h there are 11 data points which is clearly a rather small 
sample. Actually, the sample size maybe too small to draw significant conclusions, 
but nevertheless, the example is still interesting. Table H] presents simulation results, 
trying to mimic the real data scenario. We generate 11 equidistant observations on the 
interval [0,5] according to (fTTI i and (|2]i, where the errors are normal with mean zero 
and standard deviation <7 = 0.05. The NTS method seems to perform somewhat better 
in this specific scenario, although the difference with ACCEL is not dramatic. 


Table 8: Means of point estimates and actual coverage of interval estimates for the 
parameters of model ( fTTI i. where the initial value ^ is considered as known. The results 
are based on 500 simulation runs. There are 11 observations given on a uniform grid 
on [0,5], and the errors are normal with zero expectation and <7 = 0.05. The ACCEL 
point estimates are given by (fTSl) : the interval estimates are defined in ([T4l i. 


Setup 



ACCEL 

Mean Coverage 

NLS 

Mean Coverage 

n=ll 

01 

8.600e-01 

8.652e-01 

0.898 

8.600e-01 

0.944 


02 

2.079e+00 

2.082e+00 

0.924 

2.078e+00 

0.948 


03 

1.624e+00 

1.626e+00 

0.930 

1.623e+00 

0.946 


In Table |9] we see the resulting point and interval estimates based on the real data, 
using the ACCEL and NLS methods. The solutions of the system (ITtI i corresponding 
to thes e parameter estimates are displayed in Eigure|5] 

As IVan Domselaar and Hemkeil (119751) . we consider the initial va lues to be known, _ 

and use the first observation of each system state for that. We note that IVan Domselaar and Hemkei 
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Table 9: Point estima tes for the parameters of model (fTTl) based on the real data ap¬ 
pears in Varah ( 1982h . We consider the initial values to be known. The ACCEL point 
estimates are given by (fT2l i: the confidence intervals were generated according to (fT4l i. 
The left and right interval points are denoted by CI(L) and CI(R), respectively. 



Point 

ACCEL 

CI(L) 

CI(R) 

Point 

NLS 

CI(L) 

CI(R) 

01 

5.139e-01 

1.607e-01 

8.672e-01 

8.609e-01 

7.508e-01 

9.709e-01 

02 

1.552eH-00 

7.877e-01 

2.317eH-00 

2.079eH-00 

1.898eH-00 

2.260eH-00 

03 

1.316eH-00 

6.066e-01 

2.026eH-00 

1.815eH-00 

1.624eH-00 

2.006eH-00 




Figure 5: The solution to (fTTl) based on the ACCEL (solid line) and the NLS (dash 
dotted line) estimates; the observ ations are ind icated by pluses. The parameters were 
estimated using the real data from IVarahl(Il982h . 


( 1975h used an initial guess in the parameter space, whil e using the di rect integral ap¬ 
proach as we do here does not requires such a guess. IVarahl (Il982h points out that 
the data are only accurate to about 10%, and thus one should not try too hard to fit 
the data closely. Specifically, he obtained different parameter estimates, depending on 
the knots position he used for the spl ine approxim ation of the trajectories. The results 
of the NLS here are similar to those Varahl ( j_982l) report as ’VDH values’ and those 
reported in Van Domselaar and Hemkeil (ll975h .~We note that the estimates obtained 
using the ACCEL method are different from those obtained using the NLS; also, the 
ACCEL confidence intervals are wider than the NLS ones. This is unlike all the exam¬ 
ples we studied above where the ACCEL and NLS methods yielded almost identical 
results. Thus, at least in this specific example, the ACCEL method seems to be more 
conservative. However, as mentioned above, the sample size is likely to be too small 
to make any meaningful conclusions. 
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4.3 a-pinene problem 


We now consider ‘Problem 8’ of iTioa and Biedei (1991). The system is given by 


= -{9i + 92)xi{t), 

^ 4(0 = 9ixi{t), 

x'-iit) = 92Xi{t)-{93+ 94.)x3{t) + 95 x 5 ( 1 ), (18) 

■^4(0 = 93x3(1), 

4(0 = 94x3(1) - 95x5(1). 


This system characterizes a reaction that describes the thermal isomerization of a- 
pinene (jci) to dipentene (X 2 ) and alloocimene (X3), which in turn yields a — and /3 — 


pyrone ne (X 4 ) and a dimer (X 5 ). The data we use are taken from Table 2 in iBox et al 
( I 973 I) . For each state of the system, the data includes only 8 obs ervations in time. 
Thus, this is not a trivial problem to deal with, a point raised also in iTioa and Biegler 
(119911) . In Table [TOl we see the resulting point and interval estimates based on the real 
data, using the ACCEL method. We do not present the NLS method since it did not 
converge in a reasonable amount of time (using straightforward optimization as we 
did in all examples abo ve). In the last column we present the estimation result from 
Tioa and Bieglei ( 199 ih . The solution of t he system (fTSl) correspo nding to the ACCEL 


estimate is displayed in Eigure|6] Unlike iTioa and Bieglen (Il99lh . our approach does 
not require to provide an initial guess in the parameter space. 


Table 10; Point estimates for the parameters of model (fTSl ) based on the real data 
appears in Box et al. ( 1973h . We consider the initial values to be known. The ACCEL 
point estimates are given by (fT2l ): the confidence intervals were generated according to 
dirt . The left and right interval points are denoted by CI(L) and CI(R), respectively. 



Point 

CI(L) 

CI(R) 

Tioa and Biegler 119911 

01 

5.869e-05 

5.771e-05 

5.967e-05 

5.926e-05 

02 

2.830e-05 

2.740e-05 

2.920e-05 

2.963e-05 

03 

1.745e-05 

1.305e-05 

2.186e-05 

2.047e-05 

04 

2.132e-04 

1.770e-04 

2.494e-04 

2.744e-04 

05 

2.137e-05 

1.037e-05 

3.236e-05 

3.997e-05 


We conducted two simulation studies, corresponding to two different measurement 
error variances. Specifically, we generate observations according to (|2l) and (fTSl) with 
the following setup: the time grid is the same as in the real data, namely 
tj € {1230,3060,4920,7800,10680,15030,22620,36420}, resulting in a total of 8 ob¬ 
servation points. Initial values are set to the observations at the first time point, 

^ = (88.35,7.3,2.3,0.4,1.75}. The errors are normal with expectation zero and stan¬ 
dard deviations 

a = ax (44.6833,36.4111,4.9570,1.6339,12.4147} corresponding to cj,-, / = 1,..., 5. 
Here, the value a is multiplied by the mean value of each state, as calculated from the 
solutions based on the real data example. In the first study we set a = 0.02, while in 
the second we take a = 0.1. We note that the variance (7^ that corresponds to a = 0.02 
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Figure 6: The solution to (fTSl) based on the ACCEL estimates; the observations are 
indicated by different symbols, corresponding to t he system state th ey represent. The 
parameters were estimated using the real data from iBox et al. ( 1973h . 


is the order of the variance that we observed in the real data example. For each sce¬ 
nario, we repeat the experimental setup 500 times and calculate the average of point 
estimates and actual coverage of the confidence intervals. We also provide the stan¬ 
dard error of the ACCEL estimator as calculated based on 500 simulations (STE), as 
well as the square root of the average of estimates of the asymptotic variance (ASYM). 
The results are presented in Table [TT] We see that the actual coverage is not too poor, 
but nevertheless deviates noticeably from the nominal level of 95%. Further, we see 
a considerable difference between estimates of the asymptotic variance and the actual 
finite sample variance as calculated based on 500 simulations. All these results are not 
surprising, if we recall that we have at hand only 8 observations for each system state. 


5 Summary 

Parameter estimation for ordinary differential equations is a challenging problem. Ex¬ 
isting methods (certainly those that are frequently applied in practice) mainly revolve 
around the least squares estimation theme. Although in theory the least squares estima¬ 
tor possesses good (asymptotic) properties, its evaluation is far from trivial in the cur¬ 
rent setting. Alternatively, smoothing-based methods have been shown to escape some 
of the computational challenges associated with the least squares estimation, but are 
typically less efficient statistically, meaning informally that they do not squeeze max¬ 
imal possible information from the observations. In this paper we have demonstrated 
on a number of real and simulated data examples that execution of a computationally 
simple one-step Newton-Raphson type procedure on a preliminary smoothing-based 
estimator leads to rather satisfactory estimation results, that are comparable to those in 
the least squares estimation. We remark that while it is well known that the choice of 
the smoothing parameter in nonparametric estimation is a rather delicate issue, here we 
use a very simple bandwidth selection procedure, but nevertheless the obtained results 
are good, so that our method does not require much fine tuning from the user. Further¬ 
more, our estimation method enables the researcher to obtain confidence intervals. In 
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Table 11: Means of point estimates and actual coverage of interval estimates for the 
parameters of model (fTSl l. where the initial value ^ is considered as known. The results 
are based on 500 simulation runs; see the experimental setup in the text. The ACCEL 
point estimates are given by (fT2l i: the interval estimates are defined in (fl4l i. Standard 
errors of the point estimates as calculated based on the 500 simulations (denoted by 
‘STE’). Square root of the average of the estimates of the asymptotic variance, over the 
500 simulations (denoted by ‘ASYM’). 


Setup 


True 

Mean 

Coverage 

STE 

ASYM 

a = 0.02 

01 

5.926e-05 

5.920e-05 

0.758 

6.539e-07 

3.913e-07 


02 

2.963e-05 

2.958e-05 

0.806 

5.246e-07 

3.615e-07 


03 

2.047e-05 

2.042e-05 

1.000 

5.789e-07 

1.815e-06 


04 

2.744e-04 

2.709e-04 

1.000 

7.847e-06 

2.099e-05 


05 

3.997e-05 

3.878e-05 

0.998 

2.793e-06 

6.060e-06 

a = 0.1 

01 

5.926e-05 

5.910e-05 

0.768 

3.265e-06 

3.026e-02 


02 

2.963e-05 

2.945e-05 

0.820 

2.669e-06 

2.717e-03 


03 

2.047e-05 

1.993e-05 

0.998 

2.755e-06 

9.746e-06 


04 

2.744e-04 

2.452e-04 

0.946 

8.382e-05 

1.406e-04 


05 

3.997e-05 

3.103e-05 

0.940 

2.569e-05 

9.688e-05 


particular, the empirical coverage of the confidence intervals we provide is good even 
for samples as small as n = 21 in the examples we considered. On the other hand, the 
approach discussed in this work requires that all the components of the solution to O 
are observed. Thus in future work we will apply it to a smooth and match method that 
does not require a fully observed system. 

All computations were carried out using Matlab. The code will be sent by the first 
author upon request. The local polynomial estimator in some of the examples was 
based on the implementation from Caol ( 2008 ). 
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